* Calibration sub-program for CGE Model with Balistreri labor market


** Notes on calibration procedure

* Step 1: Load in raw data

* Step 2: Balance data to be consistent with SAM

* Step 3: Calibrate CES nest parameters

parameter
io0(ind,ind2)
iod0(ind,ind2)
iof0(ind,ind2)
c0(ind)
c_d0(ind)
c_f0(ind)
l0(ind)
lg0
ltot0
y0(ind)
e0(ind)
g0(ind)
g_d0(ind)
g_f0(ind)
yg0
exports0(ind)
;

$ontext
$CALL GDXXRW.EXE c22.xlsx par=c0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE c_f22.xlsx par=c_f0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE g22.xlsx par=g0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE g_f22.xlsx par=g_f0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE l22.xlsx par=l0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE io22r.xlsx par=io0 rng=A1:W23 cdim = 1 rdim = 1
$CALL GDXXRW.EXE iof22.xlsx par=iof0 rng=A1:W23 cdim = 1 rdim = 1
$CALL GDXXRW.EXE e22.xlsx par=e0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE ex22.xlsx par=exports0 rng=A1:V2 cdim = 1
$offtext

$GDXIN io22r.gdx
$LOAD io0
$GDXIN

$GDXIN iof22.gdx
$LOAD iof0
$GDXIN

$GDXIN c22.gdx
$LOAD c0
$GDXIN

$GDXIN c_f22.gdx
$LOAD c_f0
$GDXIN

$GDXIN l22.gdx
$LOAD l0
$GDXIN

$GDXIN g22.gdx
$LOAD g0
$GDXIN

$GDXIN g_f22.gdx
$LOAD g_f0
$GDXIN

$GDXIN e22.gdx
$LOAD e0
$GDXIN

$GDXIN ex22.gdx
$LOAD exports0
$GDXIN

c_d0(ind) = c0(ind) - c_f0(ind);
iod0(ind,ind2) = io0(ind,ind2) - iof0(ind,ind2);
lg0 = 1613.3;
ltot0 = sum(ind,l0(ind))+lg0;
g_d0(ind) = g0(ind) - g_f0(ind);

y0(ind) = c_d0(ind) + g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));


** Scale consumption to get aggregate balance

parameter
ratiog
ratioc
;
ratiog = sum(ind,g0(ind))/sum(ind,c0(ind));

ratioc = sum(ind,l0(ind))/((1+ratiog)*sum(ind,c0(ind)));
c0(ind) = c0(ind)*ratioc;
c_d0(ind) = c_d0(ind)*ratioc;
c_f0(ind) = c_f0(ind)*ratioc;


g0(ind) = g0(ind)*ratioc;
g_f0(ind) = g_f0(ind)*ratioc;
g_d0(ind) = g0(ind) - g_f0(ind);

y0(ind) = c_d0(ind)+  g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
yg0 = sum(ind,g0(ind))+lg0;

display
ratiog
ratioc
g0
;





** Scale exports to get trade balance

parameter
imports0
ratiot
tbalance0
;

imports0(ind) = c_f0(ind) + g_f0(ind) + sum(ind2,iof0(ind,ind2));
ratiot = sum(ind,imports0(ind))/sum(ind,exports0(ind));
exports0(ind) = exports0(ind)*ratiot;
tbalance0 = sum(ind,exports0(ind)) - (sum(ind,c_f0(ind)) + sum(ind,g_f0(ind)) + sum(ind2,sum(ind,iof0(ind,ind2))));

display
ratiot
tbalance0
exports0
imports0
;

** Scale IO matrix to get industry balance

parameter
tfd
tfs
ex
extot
;

tfd(ind) = c_d0(ind)+  g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
tfs(ind) = l0(ind) + sum(ind2,iod0(ind2,ind)+iof0(ind2,ind));
ex(ind) = tfd(ind)-tfs(ind);
extot = sum(ind,ex(ind));

display
tfd
tfs
ex
extot
;

$INCLUDE SCALE_TRADE.gms

y0(ind) = c_d0(ind) + g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
tfd(ind) = y0(ind);
tfs(ind) = l0(ind) + sum(ind2,iod0(ind2,ind)+iof0(ind2,ind));
ex(ind) = tfd(ind)-tfs(ind);
io0(ind,ind2) = iod0(ind,ind2) + iof0(ind,ind2);

display
tfd
tfs
ex
;

* Calibrate emissions factor
parameter
f_y0 "oil/gas intensity"
etot0
;


f_y0(secondary) = sum(oilgas,io0(oilgas,secondary))/y0(secondary);

mu_e(coal,elec) = 1.350/io0(coal,elec);
mu_e(coal,industrial) = .127/sum(industrial2,io0(coal,industrial2));
mu_e(coal,commercial) = .003/sum(commercial2,io0(coal,commercial2));
mu_e(oilgas,refine) = 2.296/(io0(oilgas,refine) - f_y0(refine)*exports0(refine) + f_y0(refine)*imports0(refine));
mu_e(oilgas,refine_n) = (1.473-.524)/(sum(refine_n2,io0(oilgas,refine_n2)) -sum(elec,io0(oilgas,elec)) - sum(ngd,f_y0(ngd)*exports0(ngd)) + sum(ngd,f_y0(ngd)*imports0(ngd)));
mu_e(oilgas,elec) = .524/(io0(oilgas,elec));

e0(ind) = sum(ie,mu_e(ie,ind)*io0(ie,ind));
etot0 = sum(ind,e0(ind)) - sum(oilgas,sum(secondary,mu_e(oilgas,secondary)*f_y0(secondary)*exports0(secondary))) + sum(oilgas,sum(secondary,mu_e(oilgas,secondary)*f_y0(secondary)*imports0(secondary)));

display
f_y0
mu_e
etot0
;


* Scale data such that sum(c0) = 1;

parameter
ctot
;

ctot = sum(ind,c0(ind));
c0(ind) = c0(ind)/ctot;
c_d0(ind) = c_d0(ind)/ctot;
c_f0(ind) = c_f0(ind)/ctot;
y0(ind) = y0(ind)/ctot;
l0(ind) = l0(ind)/ctot;
lg0 = lg0/ctot;
g0(ind) = g0(ind)/ctot;
g_d0(ind) = g_d0(ind)/ctot;
g_f0(ind) = g_f0(ind)/ctot;
io0(ind,ind2) = io0(ind,ind2)/ctot;
iod0(ind,ind2) = iod0(ind,ind2)/ctot;
iof0(ind,ind2) = iof0(ind,ind2)/ctot;
e0(ind) = e0(ind)/ctot;
etot0 = etot0/ctot;
exports0(ind) = exports0(ind)/ctot;
imports0(ind) = imports0(ind)/ctot;
yg0 = yg0/ctot;

* Inpute Foreign Data

parameter
ratiof
sharec
shareg
shareio
l0_f
lg0_f
g0_f
g_d0_f
g_f0_f
io0_f
iod0_f
iof0_f
ctot_f
c0_f
c_d0_f
c_f0_f
imports0_f
exports0_f
y0_f
yg0_f
;

ratiof = 3;

sharec(ind) =  c0(ind)/(c0(ind)+g0(ind)+sum(ind2,io0(ind,ind2)));
shareg(ind) =  g0(ind)/(c0(ind)+g0(ind)+sum(ind2,io0(ind,ind2)));
shareio(ind,ind2) = io0(ind,ind2)/sum(ind3,io0(ind,ind3));


l0_f(ind) = ratiof*l0(ind);
lg0_f = ratiof*lg0;

g0_f(ind) = ratiof*g0(ind);
g_f0_f(ind) = exports0(ind)*shareg(ind);
g_d0_f(ind) = g0_f(ind) - g_f0_f(ind);

io0_f(ind,ind2) = ratiof*io0(ind,ind2);
iof0_f(ind,ind2) = exports0(ind)*(1-sharec(ind)-shareg(ind))*shareio(ind,ind2);
iod0_f(ind,ind2) = io0_f(ind,ind2) - iof0_f(ind,ind2);

ctot_f = ratiof*ctot;
c0_f(ind) = ratiof*c0(ind);
c_f0_f(ind) = exports0(ind)*sharec(ind);
c_d0_f(ind) =  c0_f(ind) - c_f0_f(ind);

imports0_f(ind) = sum(ind2,iof0_f(ind,ind2)) + c_f0_f(ind) + g_f0_f(ind);
exports0_f(ind) = sum(ind2,iof0(ind,ind2)) + c_f0(ind) + g_f0(ind);

y0_f(ind) = c_d0_f(ind) + g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
yg0_f = sum(ind,g0_f(ind)) + lg0_f;

display
exports0
imports0_f
;

** Scale IO matrix to get industry balance

parameter
tfd_f
tfs_f
ex_f
extot_f
;

tfd_f(ind) = c_d0_f(ind)+  g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
tfs_f(ind) = l0_f(ind) + sum(ind2,iod0_f(ind2,ind)+iof0_f(ind2,ind));
ex_f(ind) = tfd_f(ind)-tfs_f(ind);
extot_f = sum(ind,ex_f(ind));

display
tfd_f
tfs_f
ex_f
extot_f
;

$INCLUDE SCALE_TRADEF.gms

y0_f(ind) = c_d0_f(ind) + g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
tfd_f(ind) = y0_f(ind);
tfs_f(ind) = l0_f(ind) + sum(ind2,iod0_f(ind2,ind)+iof0_f(ind2,ind));
ex_f(ind) = tfd_f(ind)-tfs_f(ind);
io0_f(ind,ind2) = iod0_f(ind,ind2) + iof0_f(ind,ind2);
yg0_f = sum(ind,g0_f(ind)) + lg0_f;

display
tfd_f
tfs_f
ex_f
;



** Calibrate labor market variables and production outer-nest

parameters
cbar_ss
cbar_f_ss
p_ss(ind)
p_f_ss
p_cbar_ss
p_cbar_f_ss
p_c_ss
p_c_f_ss
p_ex_ss
exch_ss
lambda_ss
lambda_f_ss
in_d(ind)
in_f(ind)
gtot0
gtot0_f
;

parameters
u_ss
l_ss(ind)
lg_ss
ltot_ss
w_ss
p_l_ss
ls_ss
alphap(ind)
gammay0(ind)
f_l0
f_v0
alphapg
gammayg0
f_lg0
f_vg0
lambdag0
grev_ss
tlump_ss
prof_ss(ind)

u_f_ss
l_f_ss(ind)
lg_f_ss
ltot_f_ss
w_f_ss
p_l_f_ss
ls_f_ss
alphap_f(ind)
gammay0_f(ind)
f_l0_f
f_v0_f
alphapg_f
gammayg0_f
f_lg0_f
f_vg0_f
lambdag0_f
grev_f_ss
tlump_f_ss
prof_f_ss(ind)

;


cbar_ss = sum(ind,c0(ind));
p_ss(ind) = 1;
p_cbar_ss = 1;
p_c_ss(ind) = 1;
p_ex_ss(ind) = 1;

exch_ss = 1;
lambda_ss = (cbar_ss**(-sigma))/p_cbar_ss;
in_d(ind) = sum(ind2,io0(ind2,ind));
gtot0 = sum(ind,g0(ind));

cbar_f_ss = sum(ind,c0_f(ind));
p_f_ss(ind) = 1;
p_cbar_f_ss = 1;
p_c_f_ss(ind) = 1;

lambda_f_ss = ((cbar_f_ss/ratiof)**(-sigma))/p_cbar_f_ss;
in_f(ind) = sum(ind2,io0_f(ind2,ind));
gtot0_f = sum(ind,g0_f(ind));

* labor market variables
u_ss = 0.05;
p_l_ss = 1;
w_ss = 1/(1-u_ss);
l_ss(ind) = l0(ind)/(w_ss*(1+tau_p0));
lg_ss = lg0/(w_ss*(1+tau_p0));
ltot_ss = sum(ind,l_ss(ind))+lg_ss;
ls_ss = ltot_ss/(1-u_ss);

u_f_ss = 0.05;
p_l_f_ss = 1;
w_f_ss = 1/(1-u_f_ss);
l_f_ss(ind) = l0_f(ind)/(w_f_ss*(1+tau_p0));
lg_f_ss = lg0_f/(w_f_ss*(1+tau_p0));
ltot_f_ss = sum(ind,l_f_ss(ind))+lg_f_ss;
ls_f_ss = ltot_f_ss/(1-u_f_ss);

display
cbar_ss
cbar_f_ss
lambda_ss
lambda_f_ss
;


* outer nest variables

* Domestic Industries
alphap(ind) =  ( w_ss*(1+tau_p0)*in_d(ind)**(rho(ind)-1) ) / ( w_ss*(1+tau_p0)*in_d(ind)**(rho(ind)-1) + l_ss(ind)**(rho(ind)-1)  );
gammay0(ind) = y0(ind)/(( alphap(ind)*l_ss(ind)**rho(ind)+(1-alphap(ind))*in_d(ind)**rho(ind) )**(1/rho(ind)));

f_l0(ind) = alphap(ind)*gammay0(ind)*( alphap(ind)* (l_ss(ind))**rho(ind) + (1-alphap(ind))*in_d(ind)**rho(ind) )**(1/rho(ind)-1)*(l_ss(ind))**(rho(ind)-1);
f_v0(ind) = (1-alphap(ind))*gammay0(ind)*( alphap(ind)* (l_ss(ind))**rho(ind) + (1-alphap(ind))*in_d(ind)**rho(ind) )**(1/rho(ind)-1)*(in_d(ind))**(rho(ind)-1);

* Domestic Government

alphapg =  ( w_ss*(1+tau_p0)*gtot0**(rhog-1) ) / ( w_ss*(1+tau_p0)*gtot0**(rhog-1) + lg_ss**(rhog-1)  );
gammayg0 = yg0/(( alphapg*lg_ss**rhog+(1-alphapg)*gtot0**rhog )**(1/rhog));

f_lg0 = alphapg*gammayg0*( alphapg* (lg_ss)**rhog + (1-alphapg)*gtot0**rhog )**(1/rhog-1)*(lg_ss)**(rhog-1);
f_vg0 = (1-alphapg)*gammayg0*( alphapg* (lg_ss)**rhog + (1-alphapg)*gtot0**rhog )**(1/rhog-1)*(gtot0)**(rhog-1);
lambdag0 = 1/f_vg0;


* Foreign Industries

alphap_f(ind) =  ( w_f_ss*(1+tau_p0)*in_f(ind)**(rho(ind)-1) ) / ( w_f_ss*(1+tau_p0)*in_f(ind)**(rho(ind)-1) + l_f_ss(ind)**(rho(ind)-1) );
gammay0_f(ind) = y0_f(ind)/(( alphap_f(ind)*l_f_ss(ind)**rho(ind)+(1-alphap_f(ind))*in_f(ind)**rho(ind) )**(1/rho(ind)));

f_l0_f(ind) = alphap_f(ind)*gammay0_f(ind)*( alphap_f(ind)* (l_f_ss(ind))**rho(ind) + (1-alphap_f(ind))*in_f(ind)**rho(ind) )**(1/rho(ind)-1)*(l_f_ss(ind))**(rho(ind)-1);
f_v0_f(ind) = (1-alphap_f(ind))*gammay0_f(ind)*( alphap_f(ind)* (l_f_ss(ind))**rho(ind) + (1-alphap_f(ind))*in_f(ind)**rho(ind) )**(1/rho(ind)-1)*(in_f(ind))**(rho(ind)-1);

* Foreign Government

alphapg_f =  ( w_f_ss*(1+tau_p0)*gtot0_f**(rhog-1) ) / ( w_f_ss*(1+tau_p0)*gtot0_f**(rhog-1) + lg_f_ss**(rhog-1) );
gammayg0_f = yg0_f/(( alphapg_f*lg_f_ss**rhog+(1-alphapg_f)*gtot0_f**rhog )**(1/rhog));

f_lg0_f = alphapg_f*gammayg0_f*( alphapg_f* (lg_f_ss)**rhog + (1-alphapg_f)*gtot0_f**rhog )**(1/rhog-1)*(lg_f_ss)**(rhog-1);
f_vg0_f = (1-alphapg_f)*gammayg0_f*( alphapg_f* (lg_f_ss)**rhog + (1-alphapg_f)*gtot0_f**rhog )**(1/rhog-1)*(gtot0_f)**(rhog-1);
lambdag0_f = 1/f_vg0_f;

* Household utility
psi = (1-tau_l0)*p_l_ss*lambda_ss*ls_ss**(-chi);
psi_f = (1-tau_l0)*p_l_f_ss*lambda_f_ss*(ls_f_ss/ratiof)**(-chi);


* Government transfers and profits
grev_ss = (tau_l0+tau_p0)*ltot_ss*w_ss;
tlump_ss = (tau_l0+tau_p0)*ltot_ss*w_ss  - sum(ind,g0(ind)) - (1+tau_p0)*w_ss*lg_ss;
prof_ss(ind) = p_ss(ind)*y0(ind) - in_d(ind) - (1+tau_p0)*w_ss*l_ss(ind);

grev_f_ss = (tau_l0+tau_p0)*ltot_f_ss*w_f_ss;
tlump_f_ss = (tau_l0+tau_p0)*ltot_f_ss*w_f_ss  - sum(ind,g0_f(ind)) - (1+tau_p0)*w_f_ss*lg_f_ss;
prof_f_ss(ind) = p_f_ss(ind)*y0_f(ind) - in_f(ind) - (1+tau_p0)*w_f_ss*l_f_ss(ind);


* Verify market clearing and balance

parameter
ex_d(ind)
bc_ss
bc_ss2
bc_f_ss
bc_f_ss2
;


ex_d(ind) = y0(ind) - c_d0(ind) - g_d0(ind) - exports0(ind) - sum(ind2,iod0(ind,ind2));
bc_ss = (1-tau_l0)*w_ss*ltot_ss + tlump_ss + sum(ind,prof_ss(ind)) - p_cbar_ss*cbar_ss;
bc_ss2 = (1-tau_l0)*p_l_ss*ls_ss + tlump_ss + sum(ind,prof_ss(ind)) - p_cbar_ss*cbar_ss;

ex_f(ind) = y0_f(ind) - c_d0_f(ind) - g_d0_f(ind) - exports0_f(ind) - sum(ind2,iod0_f(ind,ind2));
bc_f_ss = (1-tau_l0)*w_f_ss*ltot_f_ss  + tlump_f_ss + sum(ind,prof_f_ss(ind)) - p_cbar_f_ss*cbar_f_ss;
bc_f_ss2 = (1-tau_l0)*p_l_f_ss*ls_f_ss  + tlump_f_ss + sum(ind,prof_f_ss(ind)) - p_cbar_f_ss*cbar_f_ss;

display
ex_d
bc_ss
bc_ss2
ex_f
bc_f_ss
bc_f_ss2
prof_ss
prof_f_ss
f_l0
f_v0
f_lg0
f_vg0
f_l0_f
f_v0_f
f_lg0_f
f_vg0_f
psi
psi_f
lambda_ss
lambda_f_ss

;



** Calibrate production inner nests

* US Domestic-foreign Composite

parameters
io_id0(ind,ind2)
io_if0(ind,ind2)
cost_df0(ind,ind2)
alphadf_d(ind,ind2)
alphadf_f(ind,ind2)
gammadf(ind,ind2)
;

loop(ind,
loop(ind2,
if(io0(ind,ind2) gt 0,

io_id0(ind,ind2) = iod0(ind,ind2)/io0(ind,ind2);
io_if0(ind,ind2) = iof0(ind,ind2)/io0(ind,ind2);
cost_df0(ind,ind2) = p_ss(ind)*io_id0(ind,ind2) + (p_f_ss(ind)/exch_ss)*io_if0(ind,ind2);


alphadf_d(ind,ind2) = (p_ss(ind)/cost_df0(ind,ind2))*(io_id0(ind,ind2))**(1/sig_df(ind));
alphadf_f(ind,ind2) = ((p_f_ss(ind)/exch_ss)/cost_df0(ind,ind2))*(io_if0(ind,ind2))**(1/sig_df(ind));

gammadf(ind,ind2) = alphadf_d(ind,ind2)+alphadf_f(ind,ind2);
alphadf_d(ind,ind2) = alphadf_d(ind,ind2)/gammadf(ind,ind2);
alphadf_f(ind,ind2) = alphadf_f(ind,ind2)/gammadf(ind,ind2);

else

io_id0(ind,ind2) = 1;
io_if0(ind,ind2) = 0;
cost_df0(ind,ind2) = 1;
alphadf_d(ind,ind2) = 1;
alphadf_f(ind,ind2) = 0;
gammadf(ind,ind2) = 1;
);
);
);

parameters
io_idg0(ind)
io_ifg0(ind)
cost_dfg0(ind)
alphadfg_d(ind)
alphadfg_f(ind)
gammadfg(ind)
;

loop(ind,
if(g0(ind) gt 0,

io_idg0(ind) = g_d0(ind)/g0(ind);
io_ifg0(ind) = g_f0(ind)/g0(ind);
cost_dfg0(ind) = p_ss(ind)*io_idg0(ind) + (p_f_ss(ind)/exch_ss)*io_ifg0(ind);


alphadfg_d(ind) = (p_ss(ind)/cost_dfg0(ind))*(io_idg0(ind))**(1/sig_df(ind));
alphadfg_f(ind) = ((p_f_ss(ind)/exch_ss)/cost_dfg0(ind))*(io_ifg0(ind))**(1/sig_df(ind));

gammadfg(ind) = alphadfg_d(ind)+alphadfg_f(ind);
alphadfg_d(ind) = alphadfg_d(ind)/gammadfg(ind);
alphadfg_f(ind) = alphadfg_f(ind)/gammadfg(ind);

else

io_idg0(ind) = 1;
io_ifg0(ind) = 0;
cost_dfg0(ind) = 1;
alphadfg_d(ind) = 1;
alphadfg_f(ind) = 0;
gammadfg(ind) = 1;
);
);


* ROW Domestic-foreign Composite

parameters
io_id0_f(ind,ind2)
io_if0_f(ind,ind2)
cost_df0_f(ind,ind2)
alphadf_d_f(ind,ind2)
alphadf_f_f(ind,ind2)
gammadf_f(ind,ind2)
;

loop(ind,
loop(ind2,
if(io0_f(ind,ind2) gt 0,

io_id0_f(ind,ind2) = iod0_f(ind,ind2)/io0_f(ind,ind2);
io_if0_f(ind,ind2) = iof0_f(ind,ind2)/io0_f(ind,ind2);
cost_df0_f(ind,ind2) = p_f_ss(ind)*io_id0_f(ind,ind2) + (p_ss(ind)*exch_ss)*io_if0_f(ind,ind2);


alphadf_d_f(ind,ind2) = (p_f_ss(ind)/cost_df0_f(ind,ind2))*(io_id0_f(ind,ind2))**(1/sig_df(ind));
alphadf_f_f(ind,ind2) = ((p_ss(ind)*exch_ss)/cost_df0_f(ind,ind2))*(io_if0_f(ind,ind2))**(1/sig_df(ind));

gammadf_f(ind,ind2) = alphadf_d_f(ind,ind2)+alphadf_f_f(ind,ind2);
alphadf_d_f(ind,ind2) = alphadf_d_f(ind,ind2)/gammadf_f(ind,ind2);
alphadf_f_f(ind,ind2) = alphadf_f_f(ind,ind2)/gammadf_f(ind,ind2);

else

io_id0_f(ind,ind2) = 1;
io_if0_f(ind,ind2) = 0;
cost_df0_f(ind,ind2) = 1;
alphadf_d_f(ind,ind2) = 1;
alphadf_f_f(ind,ind2) = 0;
gammadf_f(ind,ind2) = 1;
);
);
);

parameters
io_idg0_f(ind)
io_ifg0_f(ind)
cost_dfg0_f(ind)
alphadfg_d_f(ind)
alphadfg_f_f(ind)
gammadfg_f(ind)
;

loop(ind,
if(g0_f(ind) gt 0,

io_idg0_f(ind) = g_d0_f(ind)/g0_f(ind);
io_ifg0_f(ind) = g_f0_f(ind)/g0_f(ind);
cost_dfg0_f(ind) = p_f_ss(ind)*io_idg0_f(ind) + (p_ss(ind)*exch_ss)*io_ifg0_f(ind);


alphadfg_d_f(ind) = (p_f_ss(ind)/cost_dfg0_f(ind))*(io_idg0_f(ind))**(1/sig_df(ind));
alphadfg_f_f(ind) = ((p_ss(ind)*exch_ss)/cost_dfg0_f(ind))*(io_ifg0_f(ind))**(1/sig_df(ind));

gammadfg_f(ind) = alphadfg_d_f(ind)+alphadfg_f_f(ind);
alphadfg_d_f(ind) = alphadfg_d_f(ind)/gammadfg_f(ind);
alphadfg_f_f(ind) = alphadfg_f_f(ind)/gammadfg_f(ind);

else

io_idg0_f(ind) = 1;
io_ifg0_f(ind) = 0;
cost_dfg0_f(ind) = 1;
alphadfg_d_f(ind) = 1;
alphadfg_f_f(ind) = 0;
gammadfg_f(ind) = 1;
);
);



* Energy Composite

parameters
io_ie0(ie,ind)
cost_e0(ind)
alphae(ie,ind)
gammae(ind)
;

io_ie0(ie,ind) = io0(ie,ind)/sum(ie2,io0(ie2,ind));
cost_e0(ind) = sum(ie,cost_df0(ie,ind)*io_ie0(ie,ind));

alphae(ie,ind) = (cost_df0(ie,ind)/cost_e0(ind))*(io_ie0(ie,ind))**(1/sig_e(ind));
gammae(ind) = sum(ie,alphae(ie,ind));
alphae(ie,ind) = alphae(ie,ind)/gammae(ind);

parameters
io_ieg0(ie)
cost_eg0
alphaeg(ie)
gammaeg
;

io_ieg0(ie) = g0(ie)/sum(ie2,g0(ie2));
cost_eg0 = sum(ie,cost_dfg0(ie)*io_ieg0(ie));

alphaeg(ie) = (cost_dfg0(ie)/cost_eg0)*(io_ieg0(ie))**(1/sig_eg);
gammaeg = sum(ie,alphaeg(ie));
alphaeg(ie) = alphaeg(ie)/gammaeg;

parameters
io_ie0_f(ie,ind)
cost_e0_f(ind)
alphae_f(ie,ind)
gammae_f(ind)
;

io_ie0_f(ie,ind) = io0_f(ie,ind)/sum(ie2,io0_f(ie2,ind));
cost_e0_f(ind) = sum(ie,cost_df0_f(ie,ind)*io_ie0_f(ie,ind));

alphae_f(ie,ind) = (cost_df0_f(ie,ind)/cost_e0_f(ind))*(io_ie0_f(ie,ind))**(1/sig_e(ind));
gammae_f(ind) = sum(ie,alphae_f(ie,ind));
alphae_f(ie,ind) = alphae_f(ie,ind)/gammae_f(ind);

parameters
io_ieg0_f(ie)
cost_eg0_f
alphaeg_f(ie)
gammaeg_f
;

io_ieg0_f(ie) = g0_f(ie)/sum(ie2,g0_f(ie2));
cost_eg0_f = sum(ie,cost_dfg0_f(ie)*io_ieg0_f(ie));

alphaeg_f(ie) = (cost_dfg0_f(ie)/cost_eg0_f)*(io_ieg0_f(ie))**(1/sig_eg);
gammaeg_f = sum(ie,alphaeg_f(ie));
alphaeg_f(ie) = alphaeg_f(ie)/gammaeg_f;



* Material Composite

parameters
io_im0(im,ind)
cost_m0(ind)
alpham(im,ind)
gammam(ind)
;

io_im0(im,ind) = io0(im,ind)/sum(im2,io0(im2,ind));
cost_m0(ind) = sum(im,cost_df0(im,ind)*io_im0(im,ind));

alpham(im,ind) = (cost_df0(im,ind)/cost_m0(ind))*(io_im0(im,ind))**(1/sig_m(ind));
gammam(ind) = sum(im,alpham(im,ind));
alpham(im,ind) = alpham(im,ind)/gammam(ind);

parameters
io_img0(im)
cost_mg0
alphamg(im)
gammamg
;

io_img0(im) = g0(im)/sum(im2,g0(im2));
cost_mg0 = sum(im,cost_dfg0(im)*io_img0(im));

alphamg(im) = (cost_dfg0(im)/cost_mg0)*(io_img0(im))**(1/sig_mg);
gammamg = sum(im,alphamg(im));
alphamg(im) = alphamg(im)/gammamg;

parameters
io_im0_f(im,ind)
cost_m0_f(ind)
alpham_f(im,ind)
gammam_f(ind)
;

io_im0_f(im,ind) = io0_f(im,ind)/sum(im2,io0_f(im2,ind));
cost_m0_f(ind) = sum(im,cost_df0_f(im,ind)*io_im0_f(im,ind));

alpham_f(im,ind) = (cost_df0_f(im,ind)/cost_m0_f(ind))*(io_im0_f(im,ind))**(1/sig_m(ind));
gammam_f(ind) = sum(im,alpham_f(im,ind));
alpham_f(im,ind) = alpham_f(im,ind)/gammam_f(ind);

parameters
io_img0_f(im)
cost_mg0_f
alphamg_f(im)
gammamg_f
;

io_img0_f(im) = g0_f(im)/sum(im2,g0_f(im2));
cost_mg0_f = sum(im,cost_dfg0_f(im)*io_img0_f(im));

alphamg_f(im) = (cost_dfg0_f(im)/cost_mg0_f)*(io_img0_f(im))**(1/sig_mg);
gammamg_f = sum(im,alphamg_f(im));
alphamg_f(im) = alphamg_f(im)/gammamg_f;



* Intermediate Input Composite

parameters
vv0(ind)
io_ev0(ind)
io_mv0(ind)
cost_v0(ind)
alphav_e(ind)
alphav_m(ind)
gammav(ind)
;


vv0(ind) = sum(ind2,io0(ind2,ind));
io_ev0(ind) = sum(ie,io0(ie,ind))/vv0(ind);
io_mv0(ind) = sum(im,io0(im,ind))/vv0(ind);
cost_v0(ind) = cost_e0(ind)*io_ev0(ind) + cost_m0(ind)*io_mv0(ind);

alphav_e(ind) = (cost_e0(ind)/cost_v0(ind))*(io_ev0(ind))**(1/sig_v(ind));
alphav_m(ind) = (cost_m0(ind)/cost_v0(ind))*(io_mv0(ind))**(1/sig_v(ind));

gammav(ind) = alphav_e(ind)+alphav_m(ind);
alphav_e(ind) = alphav_e(ind)/gammav(ind);
alphav_m(ind) = alphav_m(ind)/gammav(ind);

parameters
vvg0
io_evg0
io_mvg0
cost_vg0
alphavg_e
alphavg_m
gammavg
;


vvg0 = sum(ind2,g0(ind2));
io_evg0 = sum(ie,g0(ie))/vvg0;
io_mvg0 = sum(im,g0(im))/vvg0;
cost_vg0 = cost_eg0*io_evg0 + cost_mg0*io_mvg0;

alphavg_e = (cost_eg0/cost_vg0)*(io_evg0)**(1/sig_vg);
alphavg_m = (cost_mg0/cost_vg0)*(io_mvg0)**(1/sig_vg);

gammavg = alphavg_e+alphavg_m;
alphavg_e = alphavg_e/gammavg;
alphavg_m = alphavg_m/gammavg;

parameters
vv0_f(ind)
io_ev0_f(ind)
io_mv0_f(ind)
cost_v0_f(ind)
alphav_e_f(ind)
alphav_m_f(ind)
gammav_f(ind)
;


vv0_f(ind) = sum(ind2,io0_f(ind2,ind));
io_ev0_f(ind) = sum(ie,io0_f(ie,ind))/vv0_f(ind);
io_mv0_f(ind) = sum(im,io0_f(im,ind))/vv0_f(ind);
cost_v0_f(ind) = cost_e0_f(ind)*io_ev0_f(ind) + cost_m0_f(ind)*io_mv0_f(ind);

alphav_e_f(ind) = (cost_e0_f(ind)/cost_v0_f(ind))*(io_ev0_f(ind))**(1/sig_v(ind));
alphav_m_f(ind) = (cost_m0_f(ind)/cost_v0_f(ind))*(io_mv0_f(ind))**(1/sig_v(ind));

gammav_f(ind) = alphav_e_f(ind)+alphav_m_f(ind);
alphav_e_f(ind) = alphav_e_f(ind)/gammav_f(ind);
alphav_m_f(ind) = alphav_m_f(ind)/gammav_f(ind);

parameters
vvg0_f
io_evg0_f
io_mvg0_f
cost_vg0_f
alphavg_e_f
alphavg_m_f
gammavg_f
;


vvg0_f = sum(ind2,g0_f(ind2));
io_evg0_f = sum(ie,g0_f(ie))/vvg0_f;
io_mvg0_f = sum(im,g0_f(im))/vvg0_f;
cost_vg0_f = cost_eg0_f*io_evg0_f + cost_mg0_f*io_mvg0_f;

alphavg_e_f = (cost_eg0_f/cost_vg0_f)*(io_evg0_f)**(1/sig_vg);
alphavg_m_f = (cost_mg0_f/cost_vg0_f)*(io_mvg0_f)**(1/sig_vg);

gammavg_f = alphavg_e_f+alphavg_m_f;
alphavg_e_f = alphavg_e_f/gammavg_f;
alphavg_m_f = alphavg_m_f/gammavg_f;



** Calibrate consumption aggregation function

* US Domestic-Foreign Consumption Good

parameter
alphac_d(ind)
alphac_f(ind)
gammacf(ind)
;



loop(ind,
if(c0(ind) gt 0,

alphac_d(ind) = (p_ss(ind)/p_c_ss(ind))*(c_d0(ind)/c0(ind))**(1/sig_c(ind));
alphac_f(ind) =  ((p_f_ss(ind)/exch_ss)/p_c_ss(ind))*(c_f0(ind)/c0(ind))**(1/sig_c(ind));

gammacf(ind) = alphac_d(ind)+alphac_f(ind);
alphac_d(ind) = alphac_d(ind)/gammacf(ind);
alphac_f(ind) = alphac_f(ind)/gammacf(ind);

else


alphac_d(ind) = 1;
alphac_f(ind) = 0;
gammacf(ind) = 1;
);
);

* ROW Domestic-Foreign Consumption Good

parameter
alphac_d_f(ind)
alphac_f_f(ind)
gammacf_f(ind)
;



loop(ind,
if(c0_f(ind) gt 0,

alphac_d_f(ind) = (p_f_ss(ind)/p_c_f_ss(ind))*(c_d0_f(ind)/c0_f(ind))**(1/sig_c_f(ind));
alphac_f_f(ind) =  ((p_ss(ind)*exch_ss)/p_c_f_ss(ind))*(c_f0_f(ind)/c0_f(ind))**(1/sig_c_f(ind));

gammacf_f(ind) = alphac_d_f(ind)+alphac_f_f(ind);
alphac_d_f(ind) = alphac_d_f(ind)/gammacf_f(ind);
alphac_f_f(ind) = alphac_f_f(ind)/gammacf_f(ind);

else


alphac_d_f(ind) = 1;
alphac_f_f(ind) = 0;
gammacf_f(ind) = 1;
);
);

* Consumption Bundle

parameter
alphac(ind)
gammac
;

alphac(ind) = (p_c_ss(ind)/p_cbar_ss)*(c0(ind)/cbar_ss)**(1/sigmac);
gammac = sum(ind,alphac(ind));
alphac(ind) = alphac(ind)/gammac;


parameter
alphac_ff(ind)
gammac_ff
;

alphac_ff(ind) = (p_c_f_ss(ind)/p_cbar_f_ss)*(c0_f(ind)/cbar_f_ss)**(1/sigmac);
gammac_ff = sum(ind,alphac_ff(ind));
alphac_ff(ind) = alphac_ff(ind)/gammac_ff;







